Machine Learning 案例手把手:$Tennessee\;Eastman\;Process\;Simulation\;Dataset$¶
透過實作田納西-伊士曼製程（Tennessee Eastman Process, TEP）資料集，我們會帶大家了解一個機器學習專案是如何從零開始實作的，培養定義問題和解決問題的能力。在這個案例中，我們會學到：
- 資料觀察、分析與視覺化
- 如何建立機器學習模型
- 如何評量模型成效與結果視覺化
- 如何改善模型成效
內容大綱¶
1. 資料科學流程（Data Science Process）¶
2. 定義問題並了解資料來源¶
3. 撰寫程式前置作業¶
4. 資料讀取¶
5. 資料的初步觀察和前處理¶
6. 將資料建構成可以給模型的資料格式¶
7. 資料切分為訓練集和驗證集¶
8. 訓練模型：從簡單的模型開始¶
9. 評估剛剛建立好的模型¶
10. 如何改進結果¶
- 資料前處理－Feature scaling
- 模型訓練：選用 Tree-based 模型
- 模型訓練：選用 Ensemble 模型
- 資料特徵工程－Principle Component Analysis（PCA）
11. 異常檢測的特殊作法¶
- PCA 故障檢測
- Ensembling
12. 簡單的小總結¶
### Tennessee Eastman Process Simulation Dataset
- 原始論文：Downs, James J., and Ernest F. Vogel. "A plant-wide industrial process control problem." Computers & chemical engineering 17.3 (1993): 245-255.
- 模擬資料集來源：https://www.kaggle.com/averkij/tennessee-eastman-process-simulation-dataset/activity
- 田納西-伊士曼製程（Tennessee Eastman Process, TEP）是由美國 Eastman 化學公司所建立的化工模型仿真平台所模擬的製程數據，為研究該類過程的故障診斷技術提供了一個實驗平台。其產生的資料具有時間遞移性、強耦合以及非線性特徵，因此難以得到精確的數學模型，使得模型驅動方法的性能不盡理想。數據驅動法中比較有代表性的如多變量統計，例如：主成分分析、因子分析、典型相關、聚類分析…等。
### The diagram of process 整個化學反應製成主要有五個階段：反應器、冷凝器、循環壓縮機、分離器、汽提塔。主要反映過程如下：
- A（g）+ C（g）+ D（g）→ G（liq） 產物 1
- A（g）+ C（g）+ E（g）→ H（liq） 產物 2
- A（g）+ E（g）→ F（liq） 副產物
- 3D（g）→ 2F（liq） 副產物
- ### Process Faults
| Variable number | Process variable | Type |
|---|---|---|
| IDV（1） | A/C feed ratio, B composition constant（stream 4） | Step |
| IDV（2） | B composition, A/C ratio constant（stream 4） | Step |
| IDV（3） | D feed temperature（stream 2） | Step |
| IDV（4） | Reactor cooling water inlet temperature | Step |
| IDV（5） | Condenser cooling water inlet temperature | Step |
| IDV（6） | A feed loss（stream 1） | Step |
| IDV（7） | C header pressure loss-reduced availability（stream 4） | Step |
| IDV（8） | A, B, C feed composition（stream 4） | Random variation |
| IDV（9） | D feed temperature（stream 2） | Random variation |
| IDV（10） | C feed temperature（stream 4） | Random variation |
| IDV（11） | Reactor cooling water inlet temperature | Random variation |
| IDV（12） | Condenser cooling water inlet temperature | Random variation |
| IDV（13） | Reaction kinetics | Slow drift |
| IDV（14） | Reactor cooling water valve | Sticking |
| IDV（15） | Condensor cooling water valve | Sticking |
| IDV（16） | Unknown | Unknown |
| IDV（17） | Unknown | Unknown |
| IDV（18） | Unknown | Unknown |
| IDV（19） | Unknown | Unknown |
| IDV（20） | Unknown | Unknown |
Problem Define¶
- ### 利用模擬資料來預測出現哪一種操作錯誤（分類問題，OK/ Faults 1/ Fault 2/ ...）
Machine Learning task¶
載入所需模組與套件¶
【程式用法】－ import package
- import package → 匯入 package 套件
- import package as p → 匯入 package 並重新命名為 p
- import package.module1 as m1 → 將 package 底下的 module1 匯入，並重新命名為 m1
- from package import modeule1 → 從 package 將 module1 匯入（module1 存在於 package 底下）
In [1]:
import numpy as np # 矩陣操作
import pandas as pd # 處理表格類型資料
import matplotlib.pyplot as plt # 視覺化
import seaborn as sns # 進階視覺化
np.set_printoptions(edgeitems=25, linewidth=150, formatter=dict(float=lambda x: "%.3g" % x)) # 設定顯示範圍
下載資料至 Colab 空間¶
In [2]:
# TEP_FaultFree_training_100run.csv
!gdown --id 1ckTubYilJSW9q9cmp89NGpaarqQlt2NF
# TEP_Fault_training_25run.csv
!gdown --id 1ktdevLBIeFUSRwparVpyAEE7cTlKZxDg
# TEP_Fault_training_10run.csv
!gdown --id 1zwiqv6GsnCn3jbrHXZjRJhd9iW16Eeig
# TEP_Fault_testing_10run.csv
!gdown --id 1URWWeZQyO3FhS0U_nPuzMkqEb0EfX03Y
資料集為模擬運行，採樣時間均為每 3 分鐘一筆，故每小時將有 20 筆監測資料。訓練集每回運行時間為 25 小時；測試集則為 48 小時。其中蒐集異常資料時，訓練集會以正常操作開始值至第 1 小時模擬錯誤操作，因此前 20 筆為正常資料，而後 480 筆為異常資料；測試集會以正常操作開始值至第 8 小時模擬錯誤操作，因此前 160 筆為正常資料，而後 800 筆為異常資料
- 正常資料－ train: 500 筆，test: 960 筆
- 異常資料－ train: 500 筆（前 20 筆為正常資料）， test: 960 筆（前 160 筆為正常資料）
檔案說明:
- TEP_FaultFree_training_100run.csv: 模擬 100 回正常操作訓練資料集
- TEP_Fault_training_10run.csv: 模擬 10 回錯誤操作訓練資料集
- TEP_Fault_testing_10run.csv: 模擬 10 回錯誤操作測試資料集
【程式用法】－ arguments of function
- function(X, a=1, b=None) → X 為必填的參數， a 預設值為 1, b 預設值為 None
- 若未填必填參數，則會報錯；若已有預設值之參數，未填則以預設值代入計算
- 順序可以調換，但需按照參數名稱給值。例如：function(X=array, a=2) 或 function(a=2, X=array)
- 若未給參數名稱，則以 function 原定參數順序給值。例如： function(array, 2) → X=array, a=2, b=None 代入計算
In [3]:
def function(X, a=1, b=None):
print(f'X={X}, a={a}, b={b}')
In [4]:
train_normal_df = pd.read_csv('TEP_FaultFree_training_100run.csv') # 將指定名稱的 csv 檔案讀入
train_fault_df = pd.read_csv('TEP_Fault_training_10run.csv')
資料檢視¶
In [5]:
train_normal_df.head(n=5) # 資料前 n 筆
Out[5]:
In [6]:
train_normal_df.info() # 資料資訊：資料格式、筆數、欄位名稱及其形態等
In [7]:
train_fault_df.head(n=5)
Out[7]:
資料說明¶
- faultNumber: 0 (正常樣本), 1~20（操作錯誤類型）
- simulationRun: 第幾回合
- sample: 回合中的第幾個樣本
xmeas_1~xmeas_41: 製程監控參數
- Continuous process measurements
| Description | Variable number | Base case value | Units |
|---|---|---|---|
| A feed（stream 1） | xmeas_1 | 0.25052 | kscmh |
| D feed（stream 2） | xmeas_2 | 3664.0 | kgh⁻¹ |
| E feed（stream 3） | xmeas_3 | 4509.3 | kgh⁻¹ |
| A and C feed（stream 4） | xmeas_4 | 9.3477 | kscmh |
| Recycle flow（stream 8） | xmeas_5 | 26.902 | kscmh |
| Reactor feed rate（stream 6） | xmeas_6 | 42.339 | kscmh |
| Reactor pressure | xmeas_7 | 2705.0 | kPa gauge |
| Reactor level | xmeas_8 | 75 | % |
| Reactor temperature | xmeas_9 | 120.4 | °C |
| Purge rate（stream 9） | xmeas_10 | 0.33712 | kscmh |
| Product separator temperature | xmeas_11 | 80.109 | °C |
| Product separator level | xmeas_12 | 50 | % |
| Product separator pressure | xmeas_13 | 2633.7 | kPa gauge |
| Product separator underflow（stream 10） | xmeas_14 | 25.16 | m³h⁻¹ |
| Stripper level | xmeas_15 | 50 | % |
| Stripper pressure | xmeas_16 | 3102.2 | kPa gauge |
| Stripper underflow（stream 11） | xmeas_17 | 22.949 | m³h⁻¹ |
| Stripper temperature | xmeas_18 | 65.731 | °C |
| Stripper steam flow | xmeas_19 | 230.31 | kgh⁻¹ |
| Compressor work | xmeas_20 | 341.43 | kW |
| Reactor cooling water outlet temperature | xmeas_21 | 94.599 | °C |
| Separator cooling water outlet temperature | xmeas_22 | 77.297 | °C |
- Sampled process measurements
| Reactor feed analysis（stream 6） | ||||
|---|---|---|---|---|
| Component | Variable number | Base case value | Units | Sampling frequency=0.1h |
| A | xmeas_23 | 32.188 | mol% | |
| B | xmeas_24 | 8.8933 | mol% | |
| C | xmeas_25 | 26.383 | mol% | |
| D | xmeas_26 | 26.383 | mol% | |
| E | xmeas_27 | 26.383 | mol% | |
| F | xmeas_28 | 26.383 | mol% | |
| Purge gas analysis（stream 9） | ||||
| Component | Variable number | Base case value | Units | Sampling frequency=0.1h |
| A | xmeas_29 | 32.958 | mol% | |
| B | xmeas_30 | 13.823 | mol% | |
| C | xmeas_31 | 23.978 | mol% | |
| D | xmeas_32 | 1.2565 | mol% | |
| E | xmeas_33 | 18.579 | mol% | |
| F | xmeas_34 | 2.2633 | mol% | |
| G | xmeas_35 | 4.8436 | mol% | |
| H | xmeas_36 | 2.2986 | mol% | |
| Product analysis（stream 11） | ||||
| Component | Variable number | Base case value | Units | Sampling frequency=0.25h |
| D | xmeas_37 | 0.01787 | mol% | |
| E | xmeas_38 | 0.83570 | mol% | |
| F | xmeas_39 | 0.09858 | mol% | |
| G | xmeas_40 | 53.724 | mol% | |
| H | xmeas_41 | 43.828 | mol% |
- xmv_1~xmv_11: 製程操作參數
| Description | Variable number | Base case value（%） | Low limit | High limit | Units |
|---|---|---|---|---|---|
| D feed flow（stream 2） | xmv_1 | 63.053 | 0 | 5811 | kgh⁻¹ |
| E feed flow（stream 3） | xmv_2 | 53.980 | 0 | 8354 | kgh⁻¹ |
| A feed flow（stream 1） | xmv_3 | 24.644 | 0 | 1.017 | kscmh |
| A and C feed flow（stream 4） | xmv_4 | 61.302 | 0 | 15.25 | kscmh |
| Compressor recycle valve | xmv_5 | 22.210 | 0 | 100 | % |
| Purge valve（stream 9） | xmv_6 | 40.064 | 0 | 100 | % |
| Separator pot liquid flow（stream 10） | xmv_7 | 38.100 | 0 | 65.71 | m³h⁻¹ |
| Stripper liquid product flow（stream 11） | xmv_8 | 46.534 | 0 | 49.10 | m³h⁻¹ |
| Stripper steam valve | xmv_9 | 47.446 | 0 | 100 | % |
| Reactor cooling water flow | xmv_10 | 41.106 | 0 | 227.1 | m³h⁻¹ |
| Condenser cooling water flow | xmv_11 | 18.114 | 0 | 272.6 | m³h⁻¹ |
對於此資料集 X 以及 y 分別是：
- X: xmeas_1, xmeas_2, xmeas_3, ..., xmv_1, xmv_2, ...
- y: faultNumber
In [8]:
train_normal_df.describe() # 資料描述：預設針對數值型態欄位統計相關數據，包含計數、平均值、標準差、最小、第一四分位距、第二四分位距、第三四分位距、最大值等。
Out[8]:
TEP_FaultFree_training_100run.csv 僅存在 faultNumber=0（Normal），而各個監測值的值域範圍差異極大。另外，不存在著欄位中僅有一種值。
In [9]:
train_fault_df.describe() # 資料描述：預設針對數值型態欄位統計相關數據，包含計數、平均值、標準差、最小、第一四分位距、第二四分位距、第三四分位距、最大值等。
Out[9]:
TEP_Fault_training_10run.csv 中 faultNumber=1~20（Process Faults）
資料視覺化¶
取一小部分的資料來做比較：Normal data v.s. Fault data
【程式用法】- slicing
- array[start:end:step] → 取用從 start（default=0）到 end（default=-1）的元素，間隔為 step（default=1）
In [10]:
a = list(range(0, 15))
print('a:\t ', a)
print('a[:]\t ', a[:])
print('a[:2]:\t ', a[:2])
print('a[2:6]:\t ', a[2:6])
print('a[1:]:\t ', a[1:])
print('a[2:8:2]:', a[2:8:2])
【程式用法】－ pd.DataFrame
- 取用其中指定名稱的欄位 → df[column_name] 或 df.column_name
- 取用其中指定索引的筆數 → df.loc[row_name, :] 或 df.iloc[row_num, :]
- df.loc[row_name, column_name]（loc: location, 填入名稱） v.s. df.iloc[row_num, column_num]（iloc: integer-location, 填入數字）
【程式用法】 - conditional statement
輸出為 boolean （布林值）：True / False
- a == b → a 是否等於 b
- a != b → a 是否不等於 b
- a in b → a 是否是 b 的元素
- a and b → a 和 b 同時成立
- a or b → a 或 b 成立
- ...
結合 array 索引 → 只取用為 True 的索引資料
In [11]:
a = np.arange(10)
print('a:\t', a)
print('a>5:\t', a>5)
print('a[a>5]:\t', a[a>5])
In [12]:
sample_train_normal = train_normal_df[train_normal_df.simulationRun==1] # 透過條件限制 simulationRun=1 取用部分資料
sample_train_fault = train_fault_df[train_fault_df.simulationRun==1]
- FaultNumber = 1 and 2
【程式用法】－ for loop
In [13]:
for idx, i in enumerate(range(6, 12)):
print(f'{idx+1} th loop: {i}')
In [14]:
plt.figure(figsize=(15, 20)) # 設定畫布大小
for idx, i in enumerate(range(0, 5)): # 查看第 i 個特徵的分布
plt.subplot(5, 1, idx+1) # 在 5 列 1 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(sample_train_normal['sample'],
sample_train_normal.iloc[:, i+3],
label='Normal') # 繪製折線圖
for j in range(2): # 查看第 j 種 process fault
plt.plot(sample_train_fault.loc[sample_train_fault.faultNumber==j+1, 'sample'],
sample_train_fault.loc[sample_train_fault.faultNumber==j+1, f'xmeas_{i+1}'],
label=f'Fault_{j+1}') # 繪製折線圖
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.title(sample_train_normal.columns[i+3]) # 設定標題
plt.tight_layout() # 自動保持子圖之間的間距
In [15]:
plt.figure(figsize=(15, 20)) # 設定畫布大小
for idx, i in enumerate(range(0, 5)): # 查看第 i 個特徵的分布
plt.subplot(5, 1, idx+1) # 在 5 列 1 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(sample_train_normal['sample'],
sample_train_normal.iloc[:, i+3],
label='Normal') # 繪製折線圖
for j in range(2, 4):
plt.plot(sample_train_fault.loc[sample_train_fault.faultNumber==(j+1), 'sample'],
sample_train_fault.loc[sample_train_fault.faultNumber==(j+1), f'xmeas_{i+1}'],
label=f'Fault_{j+1}') # 繪製折線圖
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.title(sample_train_normal.columns[i+3]) # 設定標題
plt.tight_layout() # 自動保持子圖之間的間距
先觀察少部分的資料：對於 faultNumber = 3 and 4 與 Normal 的前 5 個特徵分布十分相近，可以再進一步觀察其他特徵的分布差異並分析此分類是否存在問題。
In [16]:
comb_df = pd.concat((sample_train_normal.iloc[20:, :], sample_train_fault[sample_train_fault['sample']>=20])) # 合併 normal 和 fault data
comb_df['faultNumber'] = comb_df['faultNumber'].astype('int') # 將欄位 faultNumber 轉換為 int（整數）型態
- 盒型圖（Boxplot）
In [17]:
plt.figure(figsize=(20, 12)) # 設定畫布大小
for idx, i in enumerate(range(3)): # 查看第 i 個特徵分布
plt.subplot(3, 1, idx+1) # 在 3 列 1 行的畫布上，選擇編號為 idx+1 的位置
sns.boxplot(data=comb_df[['faultNumber', f'xmeas_{i+1}']],
x='faultNumber',
y=f'xmeas_{i+1}') # 繪製盒鬚圖
plt.tight_layout() # 自動保持子圖之間的間距
- 小提琴圖（Violinplot）
In [18]:
plt.figure(figsize=(25, 12)) # 設定畫布大小
for idx, i in enumerate(range(3)): # 查看第 i 個特徵分布
plt.subplot(3, 1, idx+1) # 在 3 列 1 行的畫布上，選擇編號為 idx+1 的位置
sns.violinplot(data=comb_df[['faultNumber', f'xmeas_{i+1}']],
x='faultNumber',
y=f'xmeas_{i+1}') # 繪製小提琴圖
plt.tight_layout() # 自動保持子圖之間的間距
- 相關係數熱圖
In [19]:
corr = comb_df.iloc[:,[0]+[i for i in range(3, 55)]].corr(method='kendall') # 計算 faultNumber 以及 xmeas, xmv 間的相關係數
plt.figure(figsize=(30, 20)) # 設定畫布大小
sns.heatmap(corr) # 繪製熱圖（heatmap）
Out[19]:
In [20]:
comb_df['faultNumber'] = comb_df['faultNumber'].astype('str') # 將欄位 faultNumber 轉換為 str（字串）型態
new_comb_df = pd.get_dummies(comb_df, columns=['faultNumber']) # 將 faultNumber 拆成 20 個欄位（One-hot encoding：屬於該類別，給值 1 反之為 0）
In [21]:
new_comb_df.head()
Out[21]:
In [22]:
corr = new_comb_df.iloc[:,2:].corr(method='kendall') # 計算各種 faultNumber 以及 xmeas, xmv 間的相關係數
plt.figure(figsize=(30, 20)) # 繪製畫布大小
sns.heatmap(corr) # 繪製熱圖
Out[22]:
6.將資料建構成可以給模型的資料格式¶
- 所有的特徵都轉換為數字 → 模型才能做運算
- 特徵不能包含目標欄位或從目標欄位做轉換而來 → 避免模型偷看答案
- 考慮實際的流程，新一筆資料進來不具有目標值時，能蒐集到的特徵是否與訓練集一致 → 測試資料需要做的資料處理與訓練資料一致
In [23]:
train_all_10run = pd.concat((train_normal_df[train_normal_df['simulationRun']<=10], train_fault_df[train_fault_df['sample']>20])) # 合併 normal 和 fault data
In [24]:
train_all_10run.describe()
Out[24]:
In [25]:
from sklearn.model_selection import train_test_split
X_train, X_valid, y_train, y_valid = train_test_split(train_all_10run.iloc[:, 3:].values, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17) # 隨機方式
#, stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
In [26]:
print(f'The shape of X_train: {X_train.shape}\t y_train: {y_train.shape}')
print(f'The shape of X_valid: {X_train.shape}\t y_valid: {y_valid.shape}')
In [27]:
unique, counts = np.unique(train_all_10run['faultNumber'].values, return_counts=True) # 計算各種 fault 有多少筆
plt.bar(unique, counts) # 繪製直條圖
Out[27]:
In [28]:
unique, counts = np.unique(y_train, return_counts=True) # 計算各種 fault 有多少筆
plt.bar(unique, counts) # 繪製直條圖
Out[28]:
【程式用法】－ Sklearn model
- 建立模型 → model = MODEL()
- 訓練模型 → model.fit(X, y)
- 預測結果 → model.predict(X)
In [29]:
from sklearn.linear_model import LogisticRegression
In [30]:
model = LogisticRegression() # 建立模型
model.fit(X_train, y_train) # 訓練模型
Out[30]:
In [31]:
from sklearn.metrics import mean_squared_error, classification_report, confusion_matrix
- 訓練集上的結果
In [32]:
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
- 驗證集上的結果
In [33]:
# model.score(X_valid, y_valid)
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
In [34]:
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣（Confusion matrix）
10.如何改進結果¶
先看訓練集，再看驗證集！
- 訓練集結果不好（Underfitting）
- 資料處理（feature scaling, feature selection, feature engineering）
- 增加模型複雜度（選擇其他模型，調整參數）
- 訓練集結果好，驗證集結果不好（Overfitting）
- 增加資料量
- 減少雜訊資料
- 降低模型複雜度
- 增加隨機性
資料前處理 － Feature scaling¶
In [35]:
from sklearn.preprocessing import StandardScaler # (X-mean(X))/std(X)
from sklearn.preprocessing import MinMaxScaler # (X-min(X))/(max(X)-min(X))
In [36]:
def data_preprocessing(df_input, train=True, sc=None):
if train: # 若是 training data，則按照輸入資料做 feature scaling
sc = StandardScaler()
# sc = MinMaxScaler()
df = sc.fit_transform(df_input)
else: # 若是 testing data，則按照 training scale 做 feature scaling
df = sc.transform(df_input)
return df, sc
In [37]:
train_all_10run.describe()
Out[37]:
In [38]:
train_transform_data, sc = data_preprocessing(train_all_10run.iloc[:, 3:]) # feature scaling for 52 features
In [39]:
pd.DataFrame(train_transform_data, columns=train_all_10run.columns[3:]).describe()
Out[39]:
In [40]:
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
In [41]:
model = LogisticRegression() # 建立模型
model.fit(X_train, y_train) # 訓練模型
Out[41]:
- 訓練集上的結果
In [42]:
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
- 驗證集上的結果
In [43]:
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
In [44]:
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣（Confusion matrix）
模型訓練：選用 Tree-based 模型¶
In [45]:
from sklearn.tree import DecisionTreeClassifier
In [46]:
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
In [47]:
model = DecisionTreeClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
Out[47]:
- 訓練集上的結果
In [48]:
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
- 驗證集上的結果
In [49]:
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
In [50]:
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣（Confusion matrix）
Overfitting¶
模型過度擬合在訓練資料集上（綠色線）達到很高的準確率，導致在未參與訓練的驗證資料集上表現不佳
抑制 Overfitting 的方法：
- 增加資料量
- 減少雜訊資料
- 降低模型複雜度
- 增加隨機性
模型訓練：選用 Ensemble 模型¶
In [51]:
from lightgbm import LGBMClassifier
In [52]:
X_train, X_valid, y_train, y_valid = train_test_split(train_transform_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
In [53]:
model = LGBMClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
Out[53]:
- 訓練集上的結果
In [54]:
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
- 驗證集上的結果
In [55]:
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
In [56]:
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣（Confusion matrix）
資料特徵工程 － Principle Component Analysis（PCA）¶
In [57]:
from sklearn.decomposition import PCA
In [58]:
pca = PCA(n_components=0.99) # n_components: int, float or ‘mle’
pca.fit(train_transform_data) # 計算共變異數矩陣並進行特徵分解
pca_transformed_data = pca.transform(train_transform_data) # pca 降維轉換
In [59]:
pca_transformed_data.shape
Out[59]:
In [60]:
pca.explained_variance_ratio_
Out[60]:
In [61]:
pca.explained_variance_ratio_.sum()
Out[61]:
In [62]:
X_train, X_valid, y_train, y_valid = train_test_split(pca_transformed_data, # 按給定比例切分資料集
train_all_10run['faultNumber'].values,
test_size=0.25, # 切分比例
random_state=17, # 隨機方式
stratify=train_all_10run['faultNumber'].values) # 按照 faultNumber 各類比例
In [63]:
model = LGBMClassifier() # 建立模型
model.fit(X_train, y_train) # 訓練模型
Out[63]:
- 訓練集上的結果
In [64]:
# model.score(X_train, y_train)
pred_train = model.predict(X_train) # 預測結果
print(classification_report(y_train, pred_train)) # 評估結果
- 驗證集上的結果
In [65]:
pred_valid = model.predict(X_valid) # 預測結果
print(classification_report(y_valid, pred_valid)) # 評估結果
In [66]:
print(confusion_matrix(y_valid, pred_valid)) # 列出混淆矩陣（Confusion matrix）
通常在做異常檢測時，蒐集到的資料集為大量的正常資料，以及少數的異常資料。因此，若以單純的分類做法，往往會因為資料不平均而導致未能檢測出異常資料。
然而在異常檢測問題上，會利用以下的想法，來訓練模型：
- 僅訓練正常資料（學習正常資料的特徵）
- 利用正常資料的特徵轉換，對所有資料做轉換，將預期異常資料無法正常轉換或者與正常資料的特徵有較大的差異
為了更貼近實際狀況－正常資料與異常資料比例差異懸殊，將會使用模擬 100 回正常訓練資料，以及 10 回異常資料來做演示
PCA 故障檢測¶
- T$^2$ 統計量
In [67]:
train_normal_df.describe()
Out[67]:
In [68]:
sc = StandardScaler()
sc_normal = sc.fit_transform(train_normal_df.iloc[:, 3:]) # 按照 normal data scale 做 feature scaling
In [69]:
pca_normal = PCA(n_components=0.99) # n_components: int, float or ‘mle’
transformed_normal_data = pca_normal.fit_transform(sc_normal) # 計算共變異數矩陣並進行特徵分解後，做降維轉換
In [70]:
pca_normal.explained_variance_ratio_
Out[70]:
In [71]:
transformed_normal_data.shape
Out[71]:
In [72]:
T2_normal_data = (transformed_normal_data**2 / pca_normal.explained_variance_).sum(axis=1) # 計算 T^2 統計量
In [73]:
sc_fault = sc.transform(train_fault_df.iloc[:, 3:]) # 按照 normal data scale 對 fault data 做 feature scaling
transformed_fault_data = pca_normal.transform(sc_fault) # 按照 normal data 的降維方式對 fault data 做降維
In [74]:
transformed_fault_data.shape
Out[74]:
In [75]:
T2_fault_data = (transformed_fault_data**2 / pca_normal.explained_variance_).sum(axis=1) # 計算 T^2 統計量
In [76]:
plt.figure(figsize=(25, 40)) # 設定畫布大小
for idx, i in enumerate(range(20)): # 查看第 i 種 fault 的 T^2 分布
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(range(len(T2_normal_data[:500])), # 繪製折線圖
T2_normal_data[:500],
label='Normal')
plt.plot(range(len(T2_fault_data[i*500:(i+1)*500])), # 繪製折線圖
T2_fault_data[i*500:(i+1)*500],
label=f'Fault_{i+1}')
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
- Q 統計量，又名 SPE（Square Prediction Error） 統計量
In [77]:
normal_project = np.einsum('ik,kf->ifk',
transformed_normal_data,
pca_normal.components_).sum(axis=2) # 計算 normal data 的投影向量
Q_normal_data = ((normal_project-(train_normal_df.iloc[:, 3:]-train_normal_df.mean()[3:]))**2).sum(axis=1) # 計算 Q 統計量
In [78]:
fault_project = np.einsum('ik,kf->ifk',
transformed_fault_data,
pca_normal.components_).sum(axis=2) # 計算 fault data 的投影向量
Q_fault_data = ((fault_project-(train_fault_df.iloc[:, 3:]-train_fault_df.mean()[3:]))**2).sum(axis=1) # 計算 Q 統計量
In [79]:
plt.figure(figsize=(25, 40)) # 設定畫布大小
for idx, i in enumerate(range(20)): # 查看第 i 種 fault 的 Q 分布
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(range(len(Q_normal_data[:500])), # 繪製折線圖
Q_normal_data[:500],
label='Normal')
plt.plot(range(len(Q_fault_data[i*500:(i+1)*500])), # 繪製折線圖
Q_fault_data[i*500:(i+1)*500],
label=f'Fault_{i+1}')
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
Ensembling¶
- 訓練 N 個模型，利用正常訓練集中的 N-1 個特徵，去預測餘下的 1 個特徵值
- 預測異常訓練集的各個特徵
- 計算預測出來的所有特徵值，與原始值的差異（使用 MSE 等指標評估距離）
In [80]:
sc = StandardScaler()
sc_normal_df = pd.DataFrame(sc.fit_transform(train_normal_df.iloc[:, 3:]), # 按照 normal data scale 做 feature scaling
columns=train_normal_df.iloc[:, 3:].columns)
In [81]:
import tqdm
from lightgbm import LGBMRegressor
In [82]:
def train(df, cols_to_predict):
models = {} # 各個模型存放位置
pbar = tqdm.tqdm_notebook(cols_to_predict)
for col in pbar:
pbar.set_description(f'Training model for {col}')
model = LGBMRegressor(learning_rate=0.1) # 建立模型
tr_x = df.drop([col],axis=1) # 取其中 N-1 個特徵作為 X
target = df[col] # 取剩餘的 1 個特徵作為 y
model.fit(X=tr_x, y=target) # 訓練模型
models[col] = model # 將訓練好的模型存放至 models
return models
def predict(models, df, cols_to_predict):
preds = []
for col in cols_to_predict:
test_x = df.drop([col], axis=1) # 取其中 N-1 個特徵作為 X
pred = models[col].predict(test_x) # 預測剩餘的 1 個特徵
preds.append(pred)
return preds
In [83]:
sc_normal_df.describe()
Out[83]:
In [84]:
features_to_predict = sc_normal_df.columns
models = train(sc_normal_df, features_to_predict) # 透過自定義的 train 訓練模型
In [85]:
def get_mse(sample, preds):
return np.square((sample.loc[:,features_to_predict] - np.transpose(preds))**2).mean(axis=1) # 計算 mean square error
In [86]:
plt.figure(figsize=(25,40)) # 設定畫布大小
plt.title('MSE for a normal and fault samples') # 設定標題
normal_Run = np.random.randint(100)+1 # 隨機取正常資料的其中一回合
print(f'Normal simulationRun = {normal_Run}')
normal_sample = pd.DataFrame(sc.transform(train_normal_df[train_normal_df.simulationRun==normal_Run].iloc[:, 3:]),
columns=train_normal_df.columns[3:])
normal_preds = predict(models, normal_sample, features_to_predict) # 透過自定義的 predict 預測結果
fault_Run = np.random.randint(10)+1 # 隨機取異常資料的其中一回合
print(f'fault simulationRun = {fault_Run}')
for idx, i in enumerate(range(20)):
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(get_mse(normal_sample, normal_preds),
label='Normal') # 得到預測結果與原始特徵之間的 MSE 值
faulty_sample = pd.DataFrame(sc.transform(train_fault_df[(train_fault_df.simulationRun==fault_Run) & (train_fault_df.faultNumber==(i+1))].iloc[:,3:]),
columns=train_fault_df.columns[3:])
faulty_preds = predict(models, faulty_sample, features_to_predict) # 透過自定義的 predict 預測結果
plt.plot(get_mse(faulty_sample, faulty_preds),
label=f'Fault_{i+1}') # 得到預測結果與原始特徵之間的 MSE 值
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
- 不同 SimulationRun 的結果
In [87]:
plt.figure(figsize=(25,40)) # 設定畫布大小
plt.title('MSE for a normal and fault samples') # 設定標題
normal_Run = np.random.randint(100)+1 # 隨機取正常資料的其中一回合
print(f'Normal simulationRun = {normal_Run}')
normal_sample = pd.DataFrame(sc.transform(train_normal_df[train_normal_df.simulationRun==normal_Run].iloc[:, 3:]),
columns=train_normal_df.columns[3:])
normal_preds = predict(models, normal_sample, features_to_predict) # 透過自定義的 predict 預測結果
fault_Run = np.random.randint(10)+1 # 隨機取異常資料的其中一回合
print(f'fault simulationRun = {fault_Run}')
for idx, i in enumerate(range(20)):
plt.subplot(10, 2, idx+1) # 在 10 列 2 行的畫布上，選擇編號為 idx+1 的位置
plt.yscale('log') # 將 y 方向的間距取 log
plt.plot(get_mse(normal_sample, normal_preds),
label='Normal') # 得到預測結果與原始特徵之間的 MSE 值
faulty_sample = pd.DataFrame(sc.transform(train_fault_df[(train_fault_df.simulationRun==fault_Run) & (train_fault_df.faultNumber==(i+1))].iloc[:,3:]),
columns=train_fault_df.columns[3:])
faulty_preds = predict(models, faulty_sample, features_to_predict) # 透過自定義的 predict 預測結果
plt.plot(get_mse(faulty_sample, faulty_preds),
label=f'Fault_{i+1}') # 得到預測結果與原始特徵之間的 MSE 值
plt.axvline(x=20, color='r', linestyle='--') # 畫垂直線
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # 加上圖例
plt.tight_layout() # 自動保持子圖之間的間距
- 大多數的 Fault 都與 Normal 的 MSE 值有明顯的差距，如何去界定之間的 threshold，可以藉由多次的實驗訂定，或者設定分布的信心水準找到對應的統計量作為分界點。
- 在 Fault_3, Fault_9, Fault_15 幾乎分布與 Normal 無差別，也有資料顯出此三種異常出現蒐集對應有誤的情況。
此資料集也關乎時間上的關聯性，由 Fault 17 所呈現的圖形明顯看出具有週期性的變化，因此在機器學習的做法當中，可以加入往前推算一段時間的所有特徵，或者時間相關的因素在特徵當中，例如：過往時間段的平均值、最大值、最小值等等。
然而在深度學習的做法當中，有專門為時間序列所設計的 RNN 模型，也是未來可以模型方面嘗試的方向。
資料¶
- 資料處理和模型訓練是反覆檢驗的
- 資料處理應盡可能的與背景知識相呼應
- 好的且乾淨的資料可以解決大多數的問題
模型¶
- 設計特殊模型時要考慮資料輸入資料格式的類型，以及輸出與標籤相應的目標
- 模型的設計也關乎所擁有的運算資源，追求精準度或者速度
- 不管是機器學習或深度學習，單純參數的調整可以提升的預測效能有限
簡而言之，決定模型性能的常常是資料的好壞。而資料的好壞又會被資料處理的方式影響。能否找到好的資料處理方式既切合商業目標，並讓模型好學習是資料科學專案的重點。